
#ifndef __neutronics_eigen_diffusion_test_problem_
#define __neutronics_eigen_diffusion_test_problem_

#include "libmesh.h"
#include "point.h"

class Problem
{
public:
  explicit Problem (unsigned int dim, unsigned int n_materials, unsigned int e_groups, bool s_mesh=false, const std::string& fname="") ;
  
  virtual ~Problem() ;
  
	Real Keff() const ;
	
	unsigned int Egroups() const ;
  
  virtual Real Diffusion(unsigned int ig, unsigned int imat) const;
  
  virtual Real Removal(unsigned int ig, unsigned int imat) const;
  
  virtual Real NuFission(unsigned int ig, unsigned int imat) const;
  
  virtual Real Fission(unsigned int ig, unsigned int imat) const;
  
  virtual Real Scattering(unsigned src_group, unsigned int dest_group, unsigned int imat) const;
  
  virtual Real EnergyPerFission(unsigned int ig) const;
  
  virtual Real PromptFissionSpectrum(unsigned int ig) const;
  
    // boundary conditions
  virtual bool DirichletBC(const Point& pt, Real& val) const = 0 ;
  
  virtual bool NeumannBC(const Point& pt, Point& val) const = 0 ;
  
  virtual bool RobinBC(const Point& pt, Real& val1, Point& val2) const = 0 ;

  virtual Real ExactSolution(const Point& pt, unsigned int var_num=0) const ;
  
  const bool structured_mesh ;
  const std::string external_mesh_file ;
	const unsigned int dimension, n_egroups, n_materials ;
  
protected:
  
	Real kinf_exact ;

	double **xsrem, **xsfiss, **xsnufiss, **xsdiff ;	// mat * eg 
	double ***xsscatt ;	// mat * eg * eg
	double *chip ;

	double lx, ly, lz ;
};

inline
Problem::Problem(unsigned int idimension, unsigned int materials, unsigned int e_groups, bool s_mesh, const std::string& fname) 
  : dimension(idimension), n_materials(materials), n_egroups(e_groups), structured_mesh(s_mesh), external_mesh_file(fname)
{
  kinf_exact = 0.0 ;  

  chip = new double[n_materials] ;
	xsrem = new double*[n_materials] ;
	xsfiss = new double*[n_materials] ;
	xsnufiss = new double*[n_materials] ;
	xsdiff = new double*[n_materials] ;
	
	xsscatt = new double**[n_materials] ;	
	
	for (unsigned int im = 0 ; im < n_materials ; im++)
	{
		xsrem[im] = new double[n_egroups] ;
		xsfiss[im] = new double[n_egroups] ;
		xsnufiss[im] = new double[n_egroups] ;
		xsdiff[im] = new double[n_egroups] ;
		
		xsscatt[im] = new double*[n_egroups] ;		
		for (unsigned int ie = 0 ; ie < n_egroups ; ie++)
			xsscatt[im][ie] = new double[n_egroups] ;
	}
}

inline
Problem::~Problem()
{
  for (unsigned int im = 0 ; im < n_materials ; im++)
	{
		delete [] xsrem[im] ;
		delete [] xsfiss[im] ;
		delete [] xsnufiss[im] ;
		delete [] xsdiff[im] ;
		
		for (unsigned int ie = 0 ; ie < n_egroups ; ie++)
			delete [] xsscatt[im][ie] ;
    
		delete [] xsscatt[im] ;
	}
	
	delete chip ;
	delete xsrem ;
	delete xsfiss ;
	delete xsnufiss ;
	delete xsdiff ;
	delete xsscatt ;  
}

inline
Real Problem::Keff() const
{
  return kinf_exact ;
}

inline
unsigned int Problem::Egroups() const
{
  return n_egroups ;
}

inline
Real Problem::ExactSolution(const Point& pt, unsigned int var_num) const
{
	return 0 ;	// default implementation
}

#endif // #define __neutronics_eigen_diffusion_test_problem_
